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ABSTRACT 

We present detailed radiative transfer simulations of the reionization history of the 
Milky Way by metal-poor globular clusters. We identify potential metal-poor globular 
cluster candidates within the Aquarius simulation using dark matter halo velocity 
dispersions. We calculate the local ionization fields via a photon-conserving, three 
dimensional non-equilibrium chemistry code and allow the model to propagate through 
to the present day. The key feature of the model is that globular cluster formation is 
suppressed if the local gas is ionized. 

We find that our spatial treatment of the ionization field leads to drastically 
different numbers and spatial distributions when compared to models where globular 
cluster formation is simply truncated at a given redshift. We investigate a range of 
plausible ionization efficiencies to determine the effect photon-rich and photon-poor 
models have on present day globular cluster properties. We find that it is possible 
for metal-poor globular clusters to have formed via the dark matter halo formation 
channel as our secondary model (delayed formation) combined with truncation at z 
= 10 produces radial distributions statistically consistent with that of the Milky Way 
metal-poor globular clusters. 

If globular clusters do indeed form within high-redshift dark matter halos, if only 
in-part, their contributions to the reionization of the local (i.e. 2 3 h~ 3 Mpc 3 centred 
on the host galaxy) volume and mass by redshift 10 could be as high as 98% and 90%, 
respectively. In our photon poorest model, this contribution drops to 60% and 50%. 
This means globular clusters are important contributors to the reionization process 
on local scales at high-redshift until more photon-rich sources dominate the photon 
budget at later times. The surviving clusters in all models have a narrow average age 
range (mean = 13.34 Gyr, a = 0.04 Gyr) consistent with current ages estimates of the 
Milky Way metal-poor globular clusters. 

We also test a simple dynamical destruction model and estimate that ~ 60% of all 
metal-poor globular clusters formed at high redshift have since been destroyed via tidal 
interactions with the host galaxy. A final extension to our model utilises the Aquarius 
merger trees whereby suppressed globular cluster descendants can potentially become 
active sources provided they reside within a neutral region of the simulation volume 
and have no active ancestors. This addition at least doubles the number of potential 
candidates but results in extended spatial distributions. 

Key words: globular cluster - radiative transfer: methods - galaxy: formation - 
cosmology: theory - intergalactic medium 



1 INTRODUCTION 

It is now well understood that the universe emerged from the 
so-called dark ages (30 < z < 1100, Rees 1997) when light 
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from the first stars and quasars ignited and radiated large 
quantities of ionizing photons into the inter-galactic medium 
(IGM). The two primary pieces of evidence we have of this 
process occurring are a) Lya forest absorption spectra to- 
wards high-redshift quasars which show an increased opacity 
by the IGM at z = 6 ( |Fan et aL|2 006 ) and b) large-angle cos- 
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mic microwave background (CMB) polarisation anisotropics 
indicating that the median redshift of reionization was ap- 
proximately z rsj 10.4 ±1.4 ( jKomatsu et al.||2009 ). 

There are thought to have been two epochs of reioniza- 
tion: one for hydrogen (7 < z < 15, see |Loeb & Barkana 
for a review 



reionized from without (consistent with Ocvirk & Aubert 



2001 



3.5, 



David- 



and one for helium (z 
sen et alT]|1996| . Though the most likely sources of helium 
reionization are quasars, the sources of hydrogen reioniza- 
tion are far less well known. Possible sources of high redshift 



reionization include Population III stars (Bromm & Larson 
Wise||2012 and references therein), Population II 



2004 and 



stars (e.g. 



Fan et al 



Sokasian et al.|2003|), Q uasars flMadau et al.|1999| 



2001| |Fan e~ al. 2006) and_ more recently, parti- 



cle decay /annihilation (Mapelli et al. 2006|. Recen t work 



by Schaerer & Charbonnel (2011) and Conroy et al. (2011) 



suggests that initial stellar masses of today's ancient glob- 
ular clusters (hereafter GCs) could have been as much as 
8-25 times higher when they first formed, reinforcing ear- 
lier conclusions that metal-poor GCs (hereafter MPGCs) 
could have significantly contributed to the reionization of 
the IGM at high redshift. [Ricotti fe Shu"il|(|2000| were among 



the first to formally suggest that GCs could have supplied a 
large quantity of the ionizing radiation but direct attempts 
to measure their contribution were carried out by Ricotti| 
( 2002| , [Ricotti fe Ostriker] ( |2004| and |Power et al .] fl2009|) 



(though Power et al.|(|2009|) were only examining X-ray bi- 



naries within GCs as means of reionizing the IGM). Though 
insightful works, all of these studies unfortunately contain 
either (i) too many uncertain parameters (e.g. star forma- 
tion efficiencies, escape fractions, photo-ionization rates) or 
(ii) lack the resolution required to accurately resolve globu- 
lar cluster formation sites within their simulations. Despite 
these limitations, all of these studies have indicated that the 
flux from the first generation of globular clusters could have 
significantly contributed to the reionization of the Milky 
Way. 

The common approach used in many previous studies 
of GC formation (e.g. |Moore et al.|2006l|Bekki|2005| was to 
essentially force a truncation redshift (usually based on dis- 
tant quasar observations) after which no more objects could 
form. This ensured that the correct satellite numbers and 
distributions were comparable to observations. Whilst this 
'reverse engineered' type technique is useful at describing 
what type of environment is required to align theory with 
observation, it does not account for the spatial inhomogene- 
ity of the ionization field whereby sources are suppressed at 



different areas and at different times. Ciardi et al. (2003), 
Furlanetto et al.| ( |2004| ), |Barkana| ( |2QQ4| ) and |Iliev et al.] 
( 2006J all clearly demonstrated that the reionization process 
of was in fact extended in time and spatially inhomogenous 
resulting in vastly different reionization times for different 
areas of the Universe. If studies of GC formation neglect the 
inhomogeneity of reionization field, their models will result 
in inaccurate formation numbers and distributions. 

Over the years, many have studied the connection be- 
tween the reionization epoch and satellite formation. Alvarez 
et al. (2009) for example, combined a Af-body simulation 
with three-dimensional reionization calculations (1 /i -1 Gpc 
width) to determine the relationship between reionization 
history and local environment. They found that on average, 
halos with mass less than 10 13 Mq were reionized internally, 
whilst almost all halos with mass greater than 10 14 Mq were 



2011). Busha et al. (2010) combined the subhalo catalogs 



from the Via Lactea II simulation with a Gpc-scale N-body 
simulation and found that by varying the reionization time 
over the range expected for Milky Way mass halos it could 
change the number of satellite galaxies by roughly 2-3 orders 
of magnitude. 



The work of Alvarez et al. ( 2009 ) and Busha et al 



(2010) were semi-analytic studies of the reionization epoch. 
Iliev et al. (2011) however, combined a cosmological sim- 



ulation with radiative transfer and followed the formation 
of the Local Group of galaxies and nearby clusters. They 
found that for photon-poor models, the Local Group could 
have been reionized by itself (photon-rich models found that 
nearby clusters reionized the Local Group externally). Lun- 
|nan et aT] ( |2012| combined three-dimensional maps of reion- 
ization (using the semi-analytic models of |Furlanetto et aL 



2004) with the initial density field of the Aquarius simula- 
tion (Springel et al. 2008). They found that the number of 



satellites depends sensitively on the reionization model, with 
a factor of 3-4 difference for a given host halo. Most recently 



and perhaps most intriguingly, Boylan-Kolchin et al. (2011) 



have found that dissipationless dark matter simulations pre- 
dict that the majority of the most massive subhaloes of the 
Milky Way are too dense to host any bright satellites at all. 
The vast majority of these works primarily focus on the ef- 
fect of patchy reionization on the satellite (dwarf) systems of 
the Milky Way. All of them however, exclude the contribu- 
tions from primordial GCs. Whilst all of these studies show 
the consequences of a patchy reionization process on galaxy 
evolution, a model incorporating both this process and all 
of the relevant reionization sources is yet to be carried out. 



Griffen et al. (2010| (hereafter, G10) recently used cos- 



mological simulations to study the reionization of the Milky 
Way by MPGCs using the Aquarius suite of simulations. 
They concluded that with a reasonable escape fraction and 
star formation rate, the primordial MPGCs of the Milky 
Way could have ionized the Milky Way by as early as z 
~ 13. They found that the UV flux from high-z GCs had 
drastic consequences for not only the spatial and dynamical 
properties of present day (z = 0) GCs, but also for their 
overall ages as well. The main caveats of their work how- 
ever, were i) they assumed that all photons emitted from 
a potential MPGC would contribute to the reionization the 
Milky Way (much of the UV flux would have escaped into 
the outer IGM), ii) they did not model the dynamical de- 
struction of the GCs once they merged with the central halo 
and iii) did not allow for delayed star formation to take place 
whereby a previously suppressed halo could reignite if the 
recombinations were high enough and/or it became neutral 
later. 

In this paper we combine the same basic methods em- 
ployed by G10 and addresses all of the major caveats. This 
is achieved by (i) modelling the patchy reionization process 
by calculating the propagation of ionization fronts directly 
using a photon-conserving, three -dimensional ray tracin g ra- 



diative transfer code (C -Ray, |Mellema et al. 2006), (ii) 



measuring the effects of dynamical destruction on present 
day (z = 0) GC properties by including the Aquarius dataset 
with the dynamical models of Baumgardt & Makino ( 2003 ) 



and (iii) includes the effects of delayed star formation by 
combining the spatial information of halos within the Aquar- 
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ius merger trees with the state of the IGM as modelled by 
C 2 -Ray. We aim to quantitatively analyse the formation of 
primordial GCs in a Milky Way type environment and de- 
termine to what extent GCs contributed to the reionization 
of the IGM and how the their contributions affected present 
day properties of metal-poor GCs. 

This paper is organised as follows: §2 describes the 
Aquarius dark-matter simulation, how we identify globular 
cluster candidates within Aquarius and the radiative trans- 
fer (RT) code used to model MPGCs as sources of reionizing 
radiation. §3 describes the extensions to our model which in- 
clude dynamical destruction and delayed star formation. §4 
describes our results, exploring several scenarios to better 
understand the consequences of high and low escape frac- 
tions, halo concentrations, dynamical destruction and de- 
layed star formation on present day (z = 0) properties of 
our GC candidate clusters. §5 discusses the consequences of 
this work in the wider context of galaxy evolution and large 
scale reionization. §6 describes the primary conclusions of 
this paper including a discussion of the caveats of this work 
and avenues for future study. 



2 METHODOLOGY 

In this section, we describe the details of our methodology. 
We begin by summarizing the relevant details of our dataset; 
the Aquarius simulation. We then describe how we iden- 
tify globular cluster candidates using the Aquarius merger 
trees and lastly, we describe how we model GCs as ionizing 
sources. 



2.1 The Aquarius simulations 

The Aquarius suite formally consists of six different simula- 
tions of Milky Way sized galaxies, one of which is the central 
dataset used in this paper. A more detailed description of 
the simulation can be found in S08, but we review the details 
relevant to this work here. 

The Aquarius suite has same cosmological initial con- 



ditions as that of the Millennium simulation ( Springel et al 



2005). The initial size of the periodic box was 100 Mpc h~ 
Mpc set in a cosmology of Q m = 0.25, Ha = 0.75, as = 0.9, 
n s = 1, and Hubble constant H — 100 h km s _1 Mpc -1 = 
73 km s _1 Mpc -1 . We are aware that the value of as = 0.9 
is higher than the currently accepted valued of as and could 
lead to an overestimate of the ages of all candidates since all 
halos will form at earlier times. In terms of reionization, this 
means that for a fixed ionizing emissivity proportional to 
the collapsed fraction, the evolution is shifted to somewhat 



earlier times, resulting in an earlier overlap epoch (Alvarez 
et al |2006| ). Studies by |Boylan-Kolchin et aT] 1 12011) predict 
that this would have a relatively minor effect. 

The Millennium simulation was searched for halos of 
roughly Milky Way mass and without massive close neigh- 



bours in the present day (z = 0). Springel et al. (2008) also 



checked that the semi-analytic modelling applied to the tar- 
get haloes predicated them to host late-type galaxies. Other- 
wise the selection was random. These were then resimulated 
using 900 3 particles in a box of dimension 100 h~ x Mpc. Af- 
ter identifying the Lagrangian region from when each halo 
formed, the mass distribution was rerun at a much higher 



spatial and mass resolution. Although coarse particles were 
used to sample the distant regions, the resolution was such 
that the tidal field was accurately resolved at all times. 

The halos labelled, 'Aq-A' to 'Aq-F', stored snapshots 
at 128 output times equally spaced in log(a), where a = 
1/(1 + z), represents the expansion factor (between redshift 
between z = 127 and z = 0). In this paper we only focus 
on the highest resolution halo available, AqA2. Dark matter 
haloes were identified using a combination of the friends- 
of- friends algorithm and SUBFIND ( Springel et al.||2001 ). As 
described in Cole et al. (2008), haloes and their substructure 



are traced through the snapshots and linked together in a 
merger tree. The smallest halo able to be resolved in the 
simulation is of order 10 5 M©. The most pertinent details 
of the Aquarius simulation used in this study are tabulated 
in Table ^ For greater detail on the simulation, see S08. 



2.2 Identifying Primordial Globular Clusters 

Independent of whether primordial GCs are significant con- 
tributors to the reionization of the IGM, the precise mech- 
anism by which they form is still largely unknown. [Peebles 
( |1984| ) was one of the first to suggest that GCs may have 
formed within extended dark matter halos at high redshift. 
Observations of thin tidal tails by |Grillmair et aL ( 1995 ) and 
Odenkirchen et al. ( 2003 ) drastically reduced the popularity 
of this model because it was found by |Moore| ([1996), that 
tidal tails should not exist if GCs live within extended halos. 
These studies however, assume that GCs at the present pe- 
riod represent the environment from which they were born. 
There is currently no solid evidence to suggest that because 
GCs today aren't observed to contain significant traces of 
dark matter, that their more massive progenitors did not 
form within dark matter either. 

More recent work has renewed interest in the dark halo- 
GC formation scenario due to a range of numerical works 
published over the past decade (|Bromm &; Clarke] [2 002 , 
Mashchenko fe Sills|2005l|Bekki fe Yahagi|20061|Bekki et"aL] 



2007||Boley et al.||2009| |Carretta et al.||2010| ). From a dy- 
namical evolution standpoint, if GCs were in fact 8-25 times 
larger when formed than observed at present (|Schaerer fe 



Charbonnel 2011), the likelihood of a giant molecular cloud 
(/sfe ~ 0-3) of rsj 10 9 Mq collapsing without being embed- 
ded within any dark halo is at this stage, small. |Bromm| 
fe Clarke| fl2002| ) and |Mashchenko fe Sillsj fl2005| ) have also 
shown that it is possible for the extended halos of GCs to 
be tidally stripped away by the present day. More recently, 



Boley et al. (2009) carried out studies of the radial distri- 
bution of bright GCs in the Milky Way and concluded that 
it is possible that they formed in biased dark matter halos 
at high redshift. The formation of multiple stellar popula- 
tions within GCs (|Bekki fe Yahagi||2006[ |Bekki et al1|2007 



and Carretta et al.||2010 ) are also consistent with the dark 
halo-GC formation channel. 

In our work, we assume that high-z GCs did form within 
the very first collapsed minihalos and due to dynamical dis- 
ruption and violent relaxation, have since lost their dark 
matter halos. This assumption also makes it possible to 
model them as sources of reionization in current dark-matter 
only, ACDM simulations relatively straightforwardly. 
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Table 1. The basic parameters of the Aquarius simulation data used in this paper. 



Name 


[M ] 


[pc] 


N hr 




M 20 o 
[Mq] 


^200 
[kpc] 


M 50 
[Mq] 


^50 

[kpc] 


N 50 


Aq-A2 


1.370 x 10 4 


65.8 


531,570,000 


75,296,170 


1.842 x 10 12 


245.88 


2.524 x 10 12 


433.52 


184,243,536 



Notes: m p is the particle mass, e is the Plummer equivalent gravitational softening length, N^ r is the number of high resolution 
particles, and N[ r the number of low resolution particles filling the rest of the volume. M200 is the virial mass of the halo, 
defined as the mass enclosed in a sphere with mean density 200 times the critical value. r2oo gives the corresponding virial 
radius. We also give the mass and radius for a sphere of overdensity 50 times the critical density, denoted as M50 and rso- 
Note that this radius encloses a mean density 200 times the background density. Finally, A^o gives the number of simulation 
particles within rso- 



Our primary model of globular cluster formation is the 
same as that found in Section 2.2 of G10. In this work how- 
ever, we focus on metal-poor GC production only but sig- 
nificantly improve on the previous work. We focus on metal- 
poor GC formation because their relative spread of ages ( De 
Angeli et al.| [2005) places them coincidentally at the ideal 
time to be significant contributors to the reionization pro- 
cess. 

We adopt a relatively simple model to identify where 
the metal-poor GC objects would form within the Aquar- 
ius simulation based on the requirements for the collapse 
of a proto-GC gas cloud. The majority of previous studies 
adopt a similar model for GC formation, but the resolution 
of the Aquarius simulation allows us to directly measure the 
primary parameter; temperature. 



Nishi (2002) and several others found that in order for 



gas clouds to efficiently cool, they must facilitate one or a 
combination of the three mechanisms; collisional excitation 
of hydrogen and helium, radiative recombination of hydro- 
gen, and bremsstrahlung. The typical cooling function for 
a proto-GC candidate reveals an immense increase in the 
cooling rate as the temperature reaches 10 4 K, critical to 
the formation of stars. 

Assuming the gas is in quasi-static equilibrium with the 
dark matter, we can use the virial theorem to relate the ID 
velocity dispersion of the dark matter halos to temperature. 
The ratio of the 1-D internal velocity dispersion of the dark 
matter subhaloes (a v ) to the inferred virial temperature of 
the gas, T v is given by; 



kT v 
limn ' 



(i) 



where ran is the mass of a hydrogen atom and \i is the mean 
molecular weight of the gas. We adopt molecular weight of 
\i — 0.58, appropriate for a fully- ionized, primordial gas. 
Using these values, we therefore assume that whenever a 
halo's velocity dispersion increased above 11.9 kms -1 (T v ~ 
10 4 K), it is a GC candidate. The smallest halo identified 
via this method contains 1619 particles, corresponding to a 
dark matter mass of 2.17 x 10 7 M©. Throughout this work, 
'MPGCs' and 'GCs' represent the exact same objects, metal- 
poor globular clusters. 



2.3 Treatment of GCs as ionizing Sources 

The Aquarius simulation allows us to identify where and 
when candidate sources (i.e. GCs) will form and our photon- 
conserving, three dimensional, non-equilibrium chemistry 



code (C 2 -Ray, Mellema et al. 2006) will ensure the radi- 
ation from these sources is treated accurately. 



2.3.1 Radiative Transfer Calculations 

C 2 -Ray is a grid-based short characteristics ray-tracing code 
which is photon conserving and traces rays away from our 
GC sources up to each cell. Photon conservation is assured 
by adopting a finite-volume approach when calculating the 
photoionization rates and by using time-averaged optical 
depths. Although helium is not modelled directly in this 
work, C 2 -Ray assumes that where and when hydrogen is 
ionized, helium is once (but not twice) ionized because of 
the similar ionization potentials of hydrogen and helium-I. 
The evolution of the ionized fraction of hydrogen, Xi, (in- 
cluding recombinations) is governed by, 



dxi 



(1 - Xi)(Yi + U e Ci) - XiUeCVB 



(2) 



where n e is the electron density, Ti is the photo-ionization 
rate, C\ is the collisional ionization rate and as is the recom- 
bination rate. The code has been tested against 9 other rival 
ray tracing codes and was found to be on par with the best 



codes available today (Iliev et al. 2006 Iliev et al. 2009) 



For more details, see Mellema et al. (2006) and references 
therein. 

A density mesh at both 256 3 and 512 3 resolutions of the 
AqA2 halo was created for each time step using a Cloud-in- 
Cell algorithm centred on the central host maintaining a 
physical volume of 6 hT 1 Mpc. This means that the cell 
width for our high resolution run is 11 hT 1 kpc/cell, sig- 
nificantly more accurate than the 195 h~ x kpc/cell used in 
Lunnan et al. (2012). Since the Aquarius simulation con- 



sists of both high and low resolution particles in the same 
volume, composite positions of all types of particles were 
used to ensure the density was properly represented at all 
spatial scales. 



2.3.2 Mapping Sources 

Potential GC formation sites were located via the prescrip- 
tion previously discussed (see Subsection 2.2) through SQL 
queries of the publicly availably Aquarius databasqj These 
were then mapped into a grid based mesh for radiative pro- 
cessing according to the run's resolution. If more than one 



1 Aquarius database; http://www.galaxy-catalogue.dur.ac.uk 
:8080/ Aquarius/ 
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halo was identified within the same cell, those halos were 
combined and a luminosity assigned following a prescription 
based on their combined mass. For the 256 3 simulation, we 
found that only eight GC candidates formed within the same 
cell as another. When the resolution was increased to 512 3 , 
no two GCs were found to have formed within the same grid 
cell (i.e. every GC was unique within its own grid cell). 

The completeness of the identified sources is not homo- 
geneous across the C 2 -Ray box. Since the Aquarius volume 
was resimulated in a Lagrangian volume, the total volume in 
which sources can be identified is less than the total volume 
of the box used by C 2 -Ray. This incompleteness near the 
edges is yet another reason to restrict analysis to only those 
sources which have merged with the central halo. Despite 
having minor source incompleteness just near the edge of 
the C 2 -Ray volume, the density field in these outer regions 
is still accurately represented since the particles which cre- 
ated the mesh span the entire volume (high-resolution and 
low-resolution) . 



2.3.3 Calculating The Total Escape Fraction, / 7 
As highlighted by 



Schaerer & Charbonnel (2011) indicates that the ionizing 



Ricotti 



( |2002| ) and |Schaerer fe Charbom 
nel| ( |2011| , the primary problem with all GCs studies is that 
there is significant uncertainty in the emissivity of GCs in 
the early universe. Exactly how many photons are produced 
in GCs over time depends largely on three parameters; what 
fraction of the gas is converted into stars / s f e , how effec- 
tive are the stars at producing ionizing photons q, which 
reach the intergalactic medium and finally, what fraction of 
the photons manage to escape the collapsing giant molec- 
ular cloud /esc- We parameterize all of these with a single 
parameter, / 7 (/ 7 = q x / s f e x / esc ) which is equivalent to 
the number of ionizing photons/atom entering the IGM (be- 
tween time steps). Since the code is sufficiently fast and the 
number of ionizing sources is relatively small, several runs 
of varying f 1 values were able to be completed. 



Baumgardt & Kroupa (2007) found that in order to 



form a bound cluster, the star formation efficiency had to 
be one third or more. If all of the gas in a halo does not end 
up in the proto-cluster, the SFE could be considerably lower. 
The range of plausible values of the SFE adopted were, 0.15 

< /sfe < 0.7. 



Yajima et al. (2011) recently combined three dimen- 
sional radiative transfer calculations and a cosmological 
smooth particle hydrodynamic simulation to study the es- 
cape fraction, / esc , of ionizing photons of high redshift galax- 
ies. They found that the escape fraction drastically increased 
towards decreasing mass with halos of mass <10 9 Mq hav- 
ing an average / eS c ~ 0.4. As described in |Ricotti (2002), 



dust extinction and large levels of radiation are absorbed 
by the molecular cloud could drastically reduce this value. 
Our adopted escape fractions, ranging from photon-poorest 
to photon-richest scenario, were 0.1 < / eS c < 0.7. 

As in G 10, we use the the Population II efficiency curve 



of Figure 2 from Tumlinson et al. ( 2004 ) (calculated using 
the STARBURST99 code of |Leitherer et al.||1999| ) to cal- 
culate the number of ionizing photons emitted per baryon. 
Averaging over a Salpeter IMF (see Appendix of Griffen 



photon output is quite independent of the IMF slope. 

Since there are relatively few observational constraints 
on the escape fraction and star formation efficiency at high 
redshift, the total escape fraction, / 7 , is largely a free pa- 
rameter with significant uncertainty. However, by spanning 
a wide range of photon efficiencies, the consequences of hav- 
ing photon rich or photon poor metal-poor globular clusters 
on the reionization of the IGM can be examined without 
having to know what each of these quantities are exactly. 
The upper and lower limits of the total number of photons 
per baryons being produced by GC candidates in this study 
can be found in Table [2] The total range of / 7 values used 
in this study, ranging from extremely photon-poor to ex- 
tremely photon-rich, are / 7 = 150, 250, 500, 700, 1000 and 
5000 photons/baryon. Our primary assumption here is that 
the baryons follow the underlying dark matter and photons 
are produced uniformly over a given time step. Any unsup- 
pressed sources will continue producing photons in the sub- 
sequent time step. Each simulation has a time step of ^20 
Myrs. 

In addition to the ionization models already described, 
we include a 'truncation model' whereby halo formation is 
suppressed beyond a fixed redshift (ztrunc)- The only other 
requirement for these halos is that they have merged with 
the central halo by the present day (z = 0). G10 originally 
used this type of suppression model and found a Zt runc = 13 
to be the most suitable of the possible truncation redshifts. 
This truncation redshift can ultimately be arbitrary as it 
simply depends on the efficiency of the sources used. We 
adopt ztrunc = 13 here for the simple purpose of illustrat- 
ing the effect inhomogeneous ionization models have on the 
numbers and distributions of GC candidates when compared 
to a truncation models. 



2.3.4 Suppression Criterion 

Throughout this work we define 'suppression' as the inabil- 
ity of candidate halos to develop into GCs because their gas 
content is ionized locally. A radiative transfer cell is consid- 
ered ionized when the ionized fraction of hydrogen is above 
a fixed threshold (xthresh). Throughout our entire study we 
set Xthresh =0.1. GC candidates residing in a cell with an 
ionized fraction above this threshold are suppressed. 

An increase in the ionized fraction threshold essentially 
decreases the strength of the suppression process. To gain a 
full understanding of the subtle differences in the strength of 
the suppression on the overall reionization process, a range 
of ionized fraction thresholds must be tested. In Appendix B 



of Iliev et al. ( 2011 ), they discuss how source suppression de- 



et al.|2010 ) gives q « 10 000. This figure could be raised by 
moving towards a more top-heavy IMF but recent work by 



pends on the Jeans mass suppression threshold. They found 
both Xthresh = 0.1 and 0.5 yield essentially the same IGM 
evolution, apart from a slight offset in time. Only for Xthresh 
= 0.9 (weak suppression) does the evolution of the IGM be- 
gin to behave very differently. However, they were studying 
the reionization process in large scale galaxy cluster environ- 
ments and so these inferences do not necessarily translate di- 
rectly to a smaller Milky Way type environment like the one 
studied here. The only way to determine the effect of higher 
ionized fraction thresholds is to do a direct recalculation of 
one of the runs presented here but with higher threshold 
values. We defer such a study to future work (see Section 
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Table 2. The photon-poorest and photon-richest scenarios con- 
strained by previous computational studies. These limits set the 
upper and lower bounds on the total number of photons/baryon 
leaving a proto-GC, / 7 , in this study. 



Parameter 


Photon-poorest 


Photon-richest 




Model 


Model 


q 


10,000 ph/b 


10,000 ph/b 


fsfe 


0.15 


0.7 


fesc 


0.1 


0.7 


fi 


150 ph/b 


-5000 ph/b 



5 for further discussion). In any case, in this work, we only 
require that a halo is for the first time above a > 11.9 km/s 
and that it resides within a cell with an ionized fraction, 
x < 0.1. 

In terms of computation time, the highest resolution 
simulations (512 3 , Model 2: delayed) were carried out in par- 
allel on 8 nodes, each using 8 cores, requiring 1.2 GB of mem- 
ory (for all sources) and 9 hours of computation time. The 
lowest resolution simulation (256 3 , Model 1: permanently 
suppressed) took 1 hour of computation time. 



3 MODEL EXTENSIONS: DELAYED STAR 
FORMATION & DYNAMICAL 
DESTRUCTION 

In this section, we detail two extensions of the previous 
model; delayed star formation and dynamical disruption. 
Whilst these processes have been studied before, it is the 
first time that the following two models have been combined 
together in the context of GC formation. 



3.1 Delayed Star Formation 

As an extension to the previous model, we made use of the 
merger trees to trace candidate GCs spatial and dynamical 
evolution throughout the AqA simulation. In many previous 
works, the suppression of candidate clusters is treated in a 
simple, straightforward manner. If a candidate satellite re- 
sides within an ionised region, it is suppressed permanently. 
We too adopt this technique for some of our models and class 
them as 'Model V or the 'permanently suppressed' model. 

On cosmological scales, such 'suppressed once-and- 
forever' type models are a good approximation since there 
are many neighbouring sources to keep most regions ionised. 
Any descendent of a suppressed halo entering one of these 
ionised regions will remain suppressed. On small scales how- 
ever, such as the one being dealt with here, there are occa- 
sions where a candidate GC, previously suppressed, may be 
able to become active again once its gas content becomes 
neutral at a later time- step. This process effectively delays 
the star formation of a GC candidate and adds previously 
suppressed candidates to the list of legitimate GCs in the 
present day (z = 0). This important process can be closely 
examined by using the merger trees of the simulation under 
study. In this work, we denote this type of model as 'Model 



2' or the 'delayed' model. These names are used interchange- 
ably though out this work. 

Our delayed model works as follows. At each snapshot, 
we add the descendent IDs of each of the active halos to a 
cumulative list. In the subsequent snapshot, candidate halos 
are checked against this list to determine if any of the po- 
tentially active sources have been active before. If a match 
is found, this halo's descendent ID is added to the list and 
it is removed from the potential new sources for that snap- 
shot. If a match is not found, the ionisation state of the cell 
that source resides in is checked and if it is neutral (or be- 
low xthresh = 0.10) it will then become active and the usual 
quantities propagated. 

In this way, Model 2 has two types of populations form- 
ing over the course of the simulation; 'first timers' (i.e. 
source which becomes active in the same snapshot they are 
identified) and 'delayed' (i.e. they become active at a later 
snapshot than the one they are identified). This also self- 
consistently excludes any halos which have any active an- 
cestors. At early times, one would expect there to be direct 
overlap of the total number of sources using either Model 1 
and Model 2 but once delayed halos add substantial ionising 
photons to the IGM at later times, this will alter the number 
and present day properties of the candidates across the two 
models. 

In summary, Model 1 (or Ml) represents a 'permen- 
tantly suppressed' scenario whereby ionization suppresses a 
GC candidate permentantly and never ignites at any subse- 
quent time step. Model 2 (or 'delayed model') incorporates 
all of the same identification methods used in Model 1 ex- 
cept that now the descendents of previously suppressed halos 
can become active if their gas content becomes sufficiently 
neutral and they have no active ancestors. 

In this work we refer to Model 1 with f 1 — 500 (where 
/esc = 0.2 and / s f e = 0.25) as our fiducial model because 
the adopted total escape fraction is within the reasonable 
range of literature estimates. Table [3] lists the two types of 
models tested and the range of photoionisation efficiencies 
adopted in this work (e.g. our fiducial model is represented 
by Ml_512_ph500; Model 1 at 512 3 resolution with f 1 = 
500). 

3.2 Dynamical Destruction 

In the second extension to our model, we address the im- 
portant process of dynamical destruction. Although we have 
excellent resolving power to locate where candidates form, 
it is impossible to locate these halo structures in the present 
day (z = 0) as the vast majority of them will have undergone 
mergers with either the host or each other. As with G10, we 
assume the collapsed baryonic gas fraction of a GC at the 
centre of each halo will follow the same trajectory as that 
of the most bound particle of the halo. At each snapshot 
we simply identify the most bound particle corresponding 
to each GC candidate. We then follow this most bound par- 
ticle through to the present day allowing us to track the GC 
formation sites over the course of the entire simulation. 

During mergers with the central halo, the formation site 
is no longer able to be tracked by subfind and so we must 
infer the trajectory of the formation site within the host 
halo based on the most bound particle alone. Baumgar dt fc| 
Makino ( 2003 ) derived a fitting formula based on a series of 
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Table 3. The simulation suites used in this study. 'Ml' 
and 'M2' denote Model 1 (permanently suppressed 
model) and Model 2 (delayed model), respectively. Our 
fiducial model is Ml_512_ph500. A lower resolution run 
(256 3 ) of our fiducial model model is also carried out 
to test for convergence to ensure our results are not 
resolution dependent (see Sub-section |5.3| . 



Name 


# cells 


f 7 [ph/b] 


Ml_512_phl50 


512 3 


150 


Ml_512_ph250 


512 3 


250 


Ml_256_ph500 


256 3 


500 


Ml_512_ph500 


512 3 


500 


Ml_512_ph700 


512 3 


700 


Ml_512_phl000 


512 3 


1000 


Ml_512_ph5000 


512 3 


5000 


M2_512_phl50 


512 3 


150 


M2_512_ph250 


512 3 


250 


M2_512_ph500 


512 3 


500 


M2_512_ph700 


512 3 


700 


M2_512_phl000 


512 3 


1000 


M2_512_ph5000 


512 3 


5000 



iV-body models to determine whether a star cluster, given 
a few basic orbital characteristics, will be disrupted during 
its orbit around a host after such an encounter takes place. 
The expression they derived (Equation 10 from |Baumgardt| 
&; Makino||2003 ) for the dissolution time was, 

i 



M 


x Rg 


( V G 


0.33 M 


kpc 


^220 km s 



(1-e), (3) 



where M is the halo mass before infall, e is the orbital ec- 
centricity, Rg is the apocentric distance of the halo (kpc), 
Vb is the circular velocity (km/s) of the host galaxy and x 
and f3 are cluster scaling law parameters depending on the 
cluster concentration. From a fit to the lifetimes of the King 
(|1966| models with dim ensionless central potential Wo =5.0, 
Baumgardt fe Makino| ( [20031 found x = 0.75, = 1.91. 

For more concentrated, Wo = 7.0 King profile clusters, 
they found a small change in the dissolution time of a clus- 
ter (see Fig. 3 of Baumgardt &; Makino 2003 for details) 
and so for the sake of simplicity we restrict our dissolution 
analysis to Wo = 5.0 King profile clusters only. In order to 
determine whether a candidate GC survives an encounter 
with the central host halo, its calculated dissolution time 
must be compared to the time it merged with the central 
host. If the time since the candidate GC merged with the 
central halo exceeds its dissolution time, then the cluster is 
disrupted. 



4 RESULTS 

Here we present the results of our simulations organised as 
follows. Section 4.1 outlines the results for Model 1. Sec- 
tion 4.2 uses the same candidates identified in Section 1 but 
applies the dynamical disruption model outlined in Section 
3.2. Finally, Section 4.3 presents the results of Model 2 which 
includes previously suppressed sources which are allowed to 
turn on given the correct environmental and antecedent con- 
ditions. The objects used in all of the following analyses have 



been restricted to only those which merge with the final 
host by z = 0. Specifically, merged means that subfind can 
no longer recognise them as two separate structures in the 
present day. The total number of objects formed without 
suppression is 310. Though there are objects which retain 
their dark matter halos, we cite G10 as a preliminary anal- 
ysis of these objects and save more detailed study of these 
survivors for future work. 



4.1 Model 1: 'suppressed permentantly' 

4-1-1 Visualisation Of The ionization Field 

We identify GCs using the method described in Section |2.2| 
and show a visual representation of these candidates ionizing 
the IGM in Figures [I] and [2] Figure [I] shows the neutral gas 
density using a fixed emissivity of f 7 = 500 photons/baryon 
at z = 14.03, 12.07, 8.90, 7.03 within a box of 6 3 h~ 3 Mpc 3 . 
This shows a clear evolution of the ionizing fronts through- 
out the IGM during the formation of the proto-Milky Way 
particularly highlighting the inhomogeneous nature of the 
ionization field. Figure [2] shows six runs at differing emissiv- 
ities all captured at the same time (z = 10.11). The emis- 
sivities range from f 7 = 150 to f 7 = 5000 and clearly illus- 
trates how higher photo-ionization efficiencies ionize larger 
volumes. Interestingly the very central regions of the host 
galaxy remain relatively neutral (R < 100 kpc). An increase 
in the local recombination rate (oc p 2 ) in this region makes 
it difficult for central region to become ionized. 



4-1.2 Formation History & Numbers 

The formation history and overall numbers in each of our 
models are shown in Figure [3] GC formation without sup- 
pression is shown by the open bins outlining each of the sup- 
pression models (shown in colour). As expected, the relative 
numbers of objects at high-redshift are similar as the forma- 
tion environments are roughly the same. As the activation of 
more sources proceeds, large quantities of ionizing photons 
are injected into the IGM resulting in greater suppression 
in the higher escape fraction models than the lower escape 
fraction models. This amount of suppression taking place is 
most significant at z < 16. There are no active sources in 
any of the models below z = 8 except for one lone source 
forming at z = 5 in the photon-poor, Ml_512_phl50 model. 

Figure [4] shows the cumulative total of GCs forming 
indicating a factor of ^3 difference in total number from the 
lowest = 150 photons/baryon) to the highest (/ 7 = 5000 
photons/baryon) ionization efficiency. 

Table ^] shows the total number of objects formed and 
suppressed over the entire simulation for different photo- 
ionization efficiencies. As expected, there is a drastic depen- 
dence of the suppression rate on the ionization efficiency. 
The number of objects which form between the photon-poor 
(/ 7 = 150) and photon-rich (/ 7 = 5000) environments in our 
models differ by a factor of 2.7. 



4.1.3 Ages 

As illustrated by the open bins in Figure|3] the identification 
of GC candidates extends down to z = 4. Through the pro- 
cess of self-ionization, the potential number of GCs drasti- 
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Figure 1. Time evolution of the ionization field. Each panel represents a spatial slice (x-y projection), 460 /i _1 kpc thick, of the ionized 
and neutral gas density from simulations of Model 1 using photo-ionization efficiencies of f 1 = 500 at z = 14.03, 12.07, 8.90 and 7.03 
respectively. The box in each panel has a comoving width of 6 h~ x Mpc. 





Figure 2. Snapshots of the ionization field for the six photo-ionization efficiencies tested in this study all taken at the same redshift (z 
= 10.11). Again, each panel represents a spatial slices (x-y projection), 460 /i -1 kpc thick, of the ionized and neutral gas density from 
simulations of Model 1. Top left to top right: f 1 = 150, 250, 500 photons/baryon and bottom left to bottom right: f 1 = 700, 1000 and 
5000 photons/baryon. The box in each panel has a comoving width of 6 h- 1 Mpc. 



cally decreases due to the ionization flux from the first group 
of active sources. This process pushes back the average age 
of each of the GC populations. The mean age (calculated 
assuming a flat cosmology) for Model 1 is shown in Table [5] 
Overall, a higher the ionization efficiency will lead to, on av- 
erage, older GCs due to the wiping out of potential sources 
at low redshift. The youngest population on average is the 
population resulting from the Ml_512_phl50 model, having 
an average age of 13.30 Gyrs whilst the oldest population on 
average is the Ml_512_phl50 model, having an average age 



of 13.42 Gyrs. Whilst each of the model GC population's 
formation is extended in redshift-space, the corresponding 
age range they occupy is quite narrow (mean = 13.39 Gyr, 
a = 0.04 Gyr). The surviving GCs in all of the models have 
mean ages consistent with the Galactic MPGC ages deter- 
mined from the Advanced Camera for Surveys (ACS) survey 
carried out by |Marm-Franch et ah] ([2009 ) (13.5 ± 1.5 Gyrs) 
and other MPGC age studies (z > 5, |Fan et aT|2006 ). 
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Figure 3. Formation histograms for each of the respective ionization efficiencies. The number of objects as a function of redshift are 
binned in whole integer redshifts for each snapshot of the simulation. For ever higher ionization efficiencies, more photons escape into the 
IGM, resulting in fewer objects forming at lower redshift (i.e. only those identified by the velocity-temperature criterion from Aquarius) 
due to to an increased level of suppression. The clear histogram represents the formation of GCs without suppression included. 



Table 4. The number of active and suppressed sources of each of 
the various escape fractions and their corresponding suppression 
rates. 



Model 


# Active 


# Suppressed 


Suppression 


Name 


GCs 


GCs 


Rate [%] 


no suppression 


310 








Ml_512_phl50 


112 


198 


0.64 


Ml_512_ph250 


103 


207 


0.67 


Ml_512_ph500 


75 


235 


0.76 


Ml_512_ph700 


72 


238 


0.77 


Ml_512_phl000 


57 


253 


0.82 


Ml_512_ph5000 


42 


268 


0.86 



Table 5. Mean formation redshift and percentile ages for our 
metal-poor GCs candidates in each model (assuming a flat cos- 
mology). 



Model 


Zmean 


Percentile Age [Gyrs] 


Name 


Formation 


10% 


50% 


90% 


no suppression 


12.25 


12.82 


13.20 


13.36 


Ml_512_phl50 


12.92 


12.87 


13.23 


13.36 


Ml_512_ph250 


13.31 


12.92 


13.25 


13.38 


Ml_512_ph500 


13.56 


12.94 


13.25 


13.39 


Ml_512_ph700 


13.75 


13.01 


13.25 


13.39 


Ml_512_phl000 


13.92 


13.01 


13.27 


13.39 


Ml_512_ph5000 


14.05 


13.01 


13.27 


13.39 
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Figure 4. Cumulative formation history for each of the respec- 
tive ionization efficiencies. Solid lines represent active sources and 
dashed lines represent suppressed sources. The solid black line 
represents the formation of all potential cluster candidates iden- 
tified within the Aq-A2 simulation. 



1-1.4 



Spatial Distribution 



Figure [5] shows the z = radial distributions of the GCs 
in our simulations in both raw number and as a fraction 
of the total number. These were obtained by following the 
most-bound particle of each of the candidates once they 
had merged with the central halo through to the present 
day (see Section 3.2 for details). The most immediate ob- 
servation of the is the lack of surviving sources in the inner 
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50 kpc. Even for the most conservative reionization model 
(Ml_512_phl50), there are ^20 too few objects to reconcile 
the radial distribution of the AqA2 halo's MPGC distribu- 
tion with that of the Milky Way's. 

Another significant feature is the difference in shape of 
the radial distributions of the inhomogenous ionization mod- 
els (i.e. models using C 2 -Ray) and the truncation model 
(i.e. using z tr unc)- The truncation model's additional GC 
numbers in the central region originate from the inclusion 
of all potential sources before z = 13 which are unavail- 
able to the models using C 2 -Ray. As Figure [3] illustrates, 
there are a number of potential sources between z = 23 - 13 
which are otherwise suppressed in the inhomogeneous ion- 
ization models. These halos, having been included in the 
candidate sample in the truncation model reveal themselves 
in the central regions of the spatial distribution at z = 0. 
This highlights the first, primary result of this work which 
is that by treating the ionization field inhomogenous man- 
ner, it results in significantly different total numbers and z 
= radial distributions of potential GCs when compared to 
a model adopting arbitrary truncation. 

A Kolmogorov Smirnov (KS) test carried out comparing 
the six C 2 -Ray distributions against the Milky Way MPGCs 
([Fe/H] < -1) found all models except the Ml_512_phl50 
and Ml_512_ph250 models are inconsistent with the Milky 
Way metal-poor GC distribution at the 1% confidence level. 
It must be restated that none of the models are designed 
to replicate the distributions of the Milky Way and these 
results only show these two of the six distributions (photon- 
rich models) are statistically consistent with the Milky 
Way's MPGC system. 

It is possible that since we are only using one Aquarius 
halo that our conclusions are susceptible to scatter. Lunnan 
et al. (2012) found that the halo-halo scatter from the six 
Aquarius halos is up to a factor of 2-3. This is an impor- 
tant caveat of this work and means that in order to get a 
proper handle on scatter of this distribution, more Milky 
Way-type galaxies need to be examined. Further discussion 
of this scatter is discussed in Section \E\ 



The time scale for dynamical friction (iChandrasekhar 



1943) is inversely proportional to mass, so the GCs will be 
more affected by this process than the dark matter parti- 
cles (which are of order 100 times less massive). For GCs of 
mass a few times 10 5 M© at radii of 10 kpc the dynamical 
friction time scale is 10 12 years, so this is unlikely to affect 
our results. 



4-1.5 Contribution To ionization Of IGM 

Figure[6]shows the enclosed mass and volume ionized of each 
of the six models within a comoving box width (centred on 
the host galaxy) spanning 1 - 6 Mpc. Within the largest 
box size, the amount of volume and mass ionized is small 
compared to the amount of volume and mass ionized in a 
small box. This is because candidate formation density is 
largest in the central most dense region of the simulation. 
Since a large box width encloses significantly more baryons 
than a small box, the total volume and mass fraction ion- 
ized is lower. As the box width decreases the amount of 
volume and mass ionized increases as expected. Similarly as 
the photo- ionization efficiency (/ 7 ) increases, the fraction 
of volume and mass of the respective enclosed volume also 
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Figure 5. Radial distribution of candidate GCs in both raw num- 
ber and normalised to the total number. The thick black line rep- 
resents the Milky Way's metal-poor GC distribution. In addition 
to the six ionization models, the black dashed line represents the 
truncation model (i.e. GCs in the sample have formed before z = 
13 and have merged with the central halo by z = 0). The legend 
also shows the results from a two-sampled Kolmogorov-Smirnov 
test carried out between the Milky Way metal-poor distribution 
and each of the respective models. Each number corresponds to 
[H, p] where H = 1 if the test rejects the null hypothesis that they 
are drawn from the same distribution at the 1% significance level 
and p is the asymptotic p- value. 



increases as expected. In the extremely photon rich model 
(Ml_512_ph5000), as much as 80% of the entire simulated is 
ionized by z = 10. 



In terms of volume, 60% and 98% of the 2 6 h 



" 3 Mpc 3 

volume centred on the AqA halo is fully ionized by z = 10 
for the photon-poorest (Ml_512_phl50) and photon-richest 
(Ml_512_ph5000) models, respectively. In terms of the to- 
tal mass ionized, 50% and 90% of the total mass within the 
same volume was ionized by z = 10 for the photon-poor and 
photon-rich models, respectively. The subsequent decrease 
in the ionized mass fraction is due to the material recom- 
bining at low redshift. The ionized volume however, remains 
relatively constant beyond z = 8 since there are fewer sources 
below this time to increase the ionized volume fraction. 

The number of ionizing photons emitted by all active 
sources is plotted in Figure [8] The most efficient of our mod- 
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Figure 7. Redshift evolution of the volume (x v ) and mass (x m ) 
fraction ionized within a 2 3 h~ 3 Mpc 3 box centred on the host 
galaxy for each of the respective photo-ionization efficiencies. 
Whilst a significant quantity of the enclosed volume is ionized 
across each of the models, recombinations at low redshift decrease 
the mass fraction ionized since the number of active sources below 
z < 10 drastically reduces. In the late reionization era, overlap by 
neighbouring ionizing sources would ensure the remaining neutral 
gas becomes ionized by z = (see |Iliev et al.|2011| ). 



els peaks at 10 photons being produced over the course 
of the simulation. The active sources in our photon poorest 
model produces approximately 10 68 ' 5 photons within the to- 
tal 6 3 h~ 3 Mpc 3 volume. 

The primary result from this analysis is that if MPGCs 
did form via the mechanism tested, then at worst, MPGCs 
contributed a non-trivial fraction 40%, Figure |7| of the 
total ionization of the IGM within a distance of ^1 Mpc 
of the Milky Way 



4.2 Extensions to the Model: Dynamical 
Disruption & Delayed GC Formation 

In this subsection we discuss extensions to the model. In 
order to fully appreciate the effect dynamical disruption and 
delayed formation have on cluster properties, we present our 
results separately. 
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Figure 8. Cumulative number of ionizing photons emitted by 
MPGC candidates within Aquarius over the course of the simu- 
lation. 



4-2.1 Dynamical Disruption 

In order to better understand the survival rates of MPGCs 
once they merge into the central host halo, we adopted 
the dynamical disruption model of Baumgar dt &; Makino| 
( 2003 ) . The details of the method of how we include dynam- 
ical destruction are discussed in Sub-section |3.2[ Using this 
method, the dissolution times for all of the candidate clus- 
ters were calculated. Objects which have a dissolution time 
less than the look-back time to when the merged with the 
host are disrupted. The accuracy of this type of modelling 
has not yet been compared with higher resolution analy- 
ses of small scale hydro dynamic systems and so we present 
these results well aware of the several caveats involved (see 
Section 



5.4 



for details) and use this preliminary work as a 
launching platform for future research. 

The sample includes all candidate clusters identified via 
the temperature-velocity dispersion threshold discussed in 
Sub-section 2.2 As shown in Figure [9] just over half (60%) 
of the clusters identified, whether they are suppressed or not, 
will undergo dynamical disruption. Further analysis of dif- 
ferent King profiles indicated only a marginal (±6%) swing 
in survival probability. 

Figure [To] presents a percentage break down of what 
happens to the sources that don't survive through to the 
present day. The only two ways of removing sources in this 
model is by dynamical disruption or suppression by ion- 
ization. With the exception of the photon-richest model 
= 5000), approximately half of the candidate GCs in 
each model are lost due to dynamical disruption and just un- 
der half are lost due to suppression via photo-ionization. The 
reason the suppressed fraction doesn't increase more signif- 
icantly as the ionization efficiency increases is because the 
sample only examines the sub-population of clusters which 
are destroyed (i.e. all candidates of the initial population 
which were removed via suppression or disruption). For ex- 
ample, in the less efficient models, more clusters are going 
to be included in the sample because fewer clusters are sup- 
pressed at lower redshift. This also means however that a 
larger percentage of these clusters are susceptible to disrup- 
tion. 
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z = 13 z = 10 z = 7 





Figure 6. The fraction of the volume (x v ) and mass (x m ) ionized at three different redshifts (z = 7, 10 and 13) plotted simultaneously 
against the six different photo-ionization efficiencies and the width of the box surrounding the host halo in comoving coordinates. 
For each of the models, at a given redshift, a box was drawn around the central host galaxy (to within the nearest cell) and the volume 
and mass fraction ionized was calculated. Clearly, if MPGCs do form via the dark halo formation channel then their contributions to 
the ionization of their formation environments is substantial. Interestingly, at low box widths, the ionized mass fraction drops. This is 
because the number of baryons drastically incr eases towards the central regions and so the ionizing photons can't penetrate the inner 
regions of the halo. As discussed in Sub-section |2. 3. 2] the source sample could be incomplete at distances greater than 2 h~ x Mpc from 
the host (i.e. box width > 4 h~ 1 Mpc) due to the fact the the high resolution volume from which the candidates are identified is smaller 
than the total box volume of C 2 -Ray (the density field is still accurately represented at all scales however). 
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Figure 9. Dissolution time of candidate MPGCs which have 
merged with the central host by z = 0. An object with a dis- 
solution time less than the look-back time it merged with the 
central host are disrupted. 



The drastic decrease in the number of potential clusters 
(excluding the effect of suppression) has significant effects on 
the spatial distributions of the potential MPGC candidates. 



Figure [TT] plots the radial distribution of each of the models 
with and without dynamical disruption. Unlike the original 
model, none of the disruption models are statistically con- 
sistent Milky Way MPGCs at a 1% significance level. 

In order to gain a more complete understanding of the 
dynamical processes occurring within the host galaxy after 
a merger takes place, a full hydrodynamic, high-resolution 
simulation would be required in order to properly determine 
the survival rates of GCs. We present these broad brush 
dynamical destruction results to show in concept, that large 
quantities of GCs will be destroyed once they merge with 
the central host galaxy and that their radial distributions 
are currently irreconcilable with the present day metal-poor 
Galactic GC population. 

4.2.2 Model #2: Delayed GC Formation 

In the results thus far, we have presented the properties 
of MPGC candidates in the context of a model which sup- 
presses halos once and permanently. This approximation is 
accurate so long as there are sufficiently high numbers of ion- 
izing photons available to keep the entire volume ionized. If 
for instance a previously suppressed halo's descendent's gas 
content became neutral such that it could collapse to form a 
GC, this candidate would be excluded in the previous model. 
As detailed in Section |3.1| we implement a delayed model 
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Figure 10. Accounting for the empty or disrupted halos. Each 
stacked bar represents the percentage fraction lost due to sup- 
pression and dynamical disruption. With the exception of the 
photon-richest model (/ 7 = 5000), approximately half are consis- 
tently lost due to disruption and just under half are lost due to 
suppression in each of the respective models. 



whereby all suppressed candidates are still potential sources 
provided they a) enter a neutral region of the IGM and b) 
have no active ancestors. Utilising the merger trees of the 
AqA2 halo enabled us to implement this model in a rela- 
tively straight forward manner. 

Figure [12] compares the formation history of the MPGC 
candidates in both the permanently suppressed and delayed 
formation models in relation to their present day and for- 
mation galactocentric radii. At early times, the source posi- 
tions between the two models are in exact agreement with 
one another (red stars within black circles) since the IGM 
hasn't begun to be significantly ionized. As the ionizing flux 
from those first objects increases over time, the number of 
suppressed candidates also increases, resulting in a num- 
ber of extra potential sources which could become active 
later if the environmental requirements are satisfied. The 
first of these sources becomes noticeable at z ^ 19 whereby 
a few previously suppressed halos become active. Due to 
the no n- linear nature of the ionizing field, these first addi- 
tional sources then lead to more significant changes at later 
times. The majority of delayed candidates activate at z < 8 
within a galactocentric distance of R < 300 h~ x kpc. This 
is a result of there being a large neutral region in the centre 
due to the recombination rate being larger than the ioniza- 
tion rate. This neutral region can also be reasonably seen 
upon a visual inspections of the central region in Figure [2] 
It is clear that previously suppressed sources are merging 
into these regions and subsequently meeting the necessary 
environmental requirements to become active. 

Figure ^3] compares the z = radial distributions of 
all of the candidate MPGCs (cumulative and normalised) 
in both the permanently suppressed and delayed formation 
models. The immediate noticeable difference is the increase 
in the number of objects in the central region (R < 40 kpc) 
in the delayed model when compared to Model 1 of the same 
ionization efficiency. As indicated in Figure [l2| the majority 
of these new candidates are coming from suppressed sources 



120 




R/kpc 

Figure 11. Model 1 (permanent suppression) radial distributions 
of candidate MPGCs in both raw number (top panel) and nor- 
malised number (bottom panel) for models with (dashed lines) 
and without (solid lines) dynamical disruption. The effect of the 
dynamical model is stark. Candidate GC populations in the cen- 
tral regions are drastically reduced resulting in far more extended 
radial distributions. For reference, the upper panel has the same 
legend as the lower panel but we exclude it for clarity. 



becoming active in within 300 kpc of the host. These delayed 
halos subsequently merge with the host galaxy and populate 
galactocentric distances between and 100 kpc. Despite the 
extra GCs are small galactocentric radio, a KS-test carried 
out on all six Model 2 distributions rejects the null hypoth- 
esis that they are drawn from the same distribution as the 
Milky Way MPGCs at the 1% significance level. 

Figure [14] compares the z = radial distributions of 
all of the candidate MPGCs in the delayed formation mod- 
els (solid lines) and the same objects truncated at z = 10 
(dashed lines). Whilst the delayed model candidates are 
far too extended in their entirety, when they are trun- 
cated their spatial distributions are more closely consistent 
with the metal-poor GCs of the Milky Way (excluding the 
Ml_512_ph5000). This suggests that if metal-poor GCs did 
form via the dark halo formation channel, then either ex- 
ternal ionization (i.e. ionization from the Local Group) or 
ionization from other unaccounted for sources (e.g. Popula- 
tion II stars) must dominate the local environment from z 
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Figure 12. Top Panel: Comparison between Model 1 (perma- 
nently suppressed, solid black circles) and Model 2 (delayed for- 
mation, red stars) using their respective comoving galactocentric 
distance at formation and their redshift of formation. Bottom 
Panel: Similar to the top panel but this time shows the galacto- 
centric distance at z = and redshift of formation of each of the 
candidates. There are extra sources forming at z < 8 in Model 2 
due to previously suppressed sources becoming active after meet- 
ing the necessary environmental conditions within R < 300 h~ 1 
kpc of the host galaxy. 



= 10 onwards in order for their radial distributions to be 
consistent with that of the Milky Way. 



5 DISCUSSION 

In this Section, we discuss the more general achievements 
and limitations of our models with a particular emphasis on 
areas of improvement. 



5.1 Overview 

In our previous work ( Griffen et al.|2010| we set out to build 



a launching platform for more detailed work. This proof- 
by-concept work proved fruitful in a few important ways. 
Through a relatively simple GC formation model, we found 
that if GCs could have formed (in-part) at the centre of 
dark halo potentials, then the contributions of metal-poor 
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Figure 13. Cumulative (top panel) and normalised (bottom 
panel) radial distributions of all photo-ionization efficiencies in 
both Model 1 (permanently suppressed, dashed) and Model 2 
(delayed formation, solid). The delayed model adds as much as 
three times more potential GC candidates than the equivalent / 7 
permanently suppressed model (Model 1). 
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Figure 14. Radial distributions of all photo-ionization efficien- 
cies in Model 2 (solid) and the same objects truncated at z = 
10 (dashed). This shows that a better match to the Milky Way 
can be obtained if further sources of ionization are added at late 
redshifts. 
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GCs stemming from these dark halos to the reionization of 
the IGM was significant. Using the cumulative mass of these 
primordial GCs, we estimated their ionizing flux and came 
to the conclusion that it would have been possible for GCs to 
have ionized the entire Milky Way by z = 13 (for f 1 — 210). 
Despite the limited success of the model, there were a few 
glaring assumptions or limitations. For example, a) 100% 
of the ionizing flux coming from each GC was assumed to 
have gone into ionizing the Milky Way, b) only one value 
for the emissivity of the GCs was examined, c) the state 
of the IGM was entirely unknown at high redshifts and d) 
an arbitrary truncation redshift based on relatively naive as- 
sumptions about the nature of the IGM. This work provided 
us with the understanding required to develop a more de- 
tailed reionization model and to look at ways of improving 
the treatment of the ionization field. 

In this latest work, we improve on all of these areas 
and more. The first immediate difference is that the ion- 
ization field emanating from the sources is inhomogeneous. 
This reveals itself in both the visuals of Figure^ and the ex- 
tended formation histograms in Figure [3] Another key fea- 
ture of Figure [3] i s that unlike the formation history in G10 



and Bekki et al. (2008), MPGC formation is extended m 
time rather than abruptly ending at one specific redshift. 
This extended formation can only come about due to the 
inhomogeneous way in which the ionization field propagates 
through the bulk IGM. Our models are consistent with the 
overall picture of extended reionization. For the photon- 
poorest model (Ml_512_phl50), the last source is activated 
as late as z = 4 whereas in the most photon-rich model 
(Ml_512_ph5000), MPGCs cease forming by z = 7 (see Ta- 
ble^. 



5.2 Globular Cluster Dark Halo Formation 
Hypothesis 

Whether or not MPGCs originally formed within dark mat- 
ter halos is still an open question. In this work we do not 
necessitate that they do but rather ask what sort of esti- 
mations can one derive if the dark matter halo formation 
channel produces distributions similar to observations. We 
conclude that although some GCs may have formed in dark 
matter halos and were subsequently stripped upon merging 
with the central host galaxy, a portion of the Galactic glob- 
ular cluster population we observe today may have come 
about by other means (e.g. mergers, giant molecular clouds, 
tidal tails). In Model 1, even before including effects like 
dynamical disruption, the radial distributions are either too 
shallow or don't contain enough GCs to reproduce the Milky 
Way distribution (even if one assumes several GCs formed 
per dark matter halo). Though the delayed model (Model 
2) does increase the number of potential sources, a great 
majority of them don't populate the required locations in 
order to reconcile the simulated distribution with the Milky 
Way MPGC distribution. This however ignores the contribu- 
tion from external sources inundating the IGM with ionizing 
photons which would suppress low redshift GCs. If we adopt 
Model 2 and truncate sources at z = 10, the spatial distri- 
butions are consistent with that of the Milky Way MPGCs. 
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Figure 15. Convergence test for our fiducial model, 
Ml_512_ph500. In terms of the cumulative formation num- 
ber, the lower resolution, 256 3 mesh shows a near-complete 
agreement with the higher resolution 512 3 mesh. A KS test 
reveals that the hypothesis that they are drawn from the same 
distribution cannot be rejected at the 5% significance level. 



5.3 Convergence 

An important concern for the models is that the results are 
not dependent on the resolution of the particular simula- 
tion used. G10 has already shown that the identification of 
sources is not resolution-limited by comparing the forma- 
tion history of objects found within the higher resolution 
AqA2 halo to the lower resolution AqA3 halo. Here, we are 
interested in the dependence of the results on the resolu- 
tion used by C 2 -Ray. In conjunction with each of the 512 3 
suite of models, a lower-resolution 256 3 simulation was car- 
ried out on our fiducial model (Ml_512_ph500) to test for 
convergence. In Figure [l5| we present the formation history 
of metal-poor GCs in our fiducial model at both resolutions, 
256 3 and 512 3 . Both models are in excellent agreement with 
one another varying only by few percent in number. A KS- 
test also verifies that they are consistent with having been 
drawn from the same distribution at the 5% significance 
level. This agreement shows that the results are not signifi- 
cantly biased by the resolution of the simulations used. 



5.4 Caveats Of This Work 

Although we have tried to eliminate as many known weak- 
nesses of the models used in this study, there are still a 
number of remaining caveats which need to be addressed in 
future studies to improve upon this work. 

• Accountability of the ionization sources; whilst the ma- 
jority of the local ionizing flux will come from low mass 
dark matter halos in the early Universe as time evolves, 
each of our models diverges from an accurate account of 
the true ionizing sources. This is because field stars, non- 
dark matter dominated giant molecular clouds and other 
such objects will begin to play a dominating role at lower 
redshifts which are not yet modelled in the current work. In- 



deed recent work by Iliev et al. (2011 ); Alvarez et al. (2009) 



suggests a reionization overlap scenario for the local Milky 
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Way volume whereby reionization was kickstarted by low 
mass sources locally, but their effective flux is superseded 
by with the ionizing flux of the Local Group. Nearby proto- 
clusters (e.g. Virgo and Fornax) clusters would also have 
had an appreciable affect at later times. When and where 
this overlap occurred is critical in understanding the origin 
of present day satellite structures. 

• Unmerged GCs; Throughout this work we have limited 
our sample to include only those which have merged with 
the central halo by z = 0. This was done to not only ensure 
an accurate sample but for simplicity. Candidates which do 
not merge with the central halo by z = could be Galactic 
dwarf galaxies as they will have retained their dark matter 
halos, though this is highly speculative. Figure 10 of G10 
shows the spatial distribution of these no n- merged objects 
is consistent with those GCs found surrounding neighbour- 
ing dwarf galaxies (see |Mateo| [l998). Overall there are 47 
distinct halos containing 68 GCs at z = 0, indicating that 
there has been not only major merging with the central host, 
but also minor- merging of GCs with proto-dwarf galaxies. 
Whilst there is sufficient research to suggest GCs could have 
formed in dark matter halos and were subsequently stripped 
of their outer envelope, it remains unclear where the degen- 
eracy breaks between GC & dwarf galaxy formation pro- 
cesses. An examination of these non-merged objects needs 
to be carried out in a more detailed manner to determine 
their connection (if any) to the known dark matter domi- 
nated satellites of Milky Way. 

• GC formation efficiency; All of the work presented here 
assumes that for every dark matter halo, exactly one GC will 
form in the centre of its potential well. Given that fragmen- 
tation of large gas clouds can lead to multiple stellar pop- 
ulations within a single dark matter halo, it still remains 
unclear as to exactly how many GCs will form in a given 
dark halo. In any case, it would only shift the halo centric 
radii plots upwards and not really have any impact on the 
structural properties of the distributions. Understanding GC 
formation efficiencies better will lead to better estimates of 
their overall numbers if they do form via the dark halo for- 
mation channel. 

• Statistical limitation using one Aquarius halo; Perhaps 
one of the biggest drawbacks of this work is that it only ex- 
amines one Milky Way type galaxy. In order to get a better 
understanding of the underlying physics involved, a larger 
sample of Milky Way type galaxies would be required. The 
statistical variation in the number of satellites in the Aquar- 
ius simulation has also been shown to vary by a factor of 2 or 
more ( |Lunnan et al.|2Q12| . The overall conclusion that GCs 
contributed a non-trivial amount to the reionization of the 
local IGM would be unlikely to change given that the most 
conservative of our models reionizes more than 50% both by 
mass and volume of the local IGM by z = 10 (see Fig. 
Although the Aquarius suite consists of six such simulations, 
the logistics behind analysing the remaining five halos was 
not feasible at this time. 

• Treatment of dynamical disruption; the exact nature of 
the dynamical processes inside the central host halo after 
a halo merges is inherently difficult to treat properly. In 
our model, we trace the single most bound particle of the 
in-falling MPGCs. These trace particles could be suscepti- 
ble to violent dynamical changes within the halo resulting in 
spatial distributions that might not necessarily reflect where 



the true z = position of the baryonic content ended up. In 
an improvement to this, one could trace the 10% most bound 
particles to get a better hold on the z = position of the 
clusters under study. Whilst we observe that the disruption 
process is a critical component in determine the final z = 
numbers (e.g. as many as half of the original GC population 
to have formed at high redshift is destroyed by dynamical 
disruption), a lot more work needs to be done to more ac- 
curately determine the degree to which GCs are destroyed 
over the course of their evolution. 

• Dependence on ionization cell-fraction threshold; Since 
our simulations do not (yet) have gas and do not calcu- 
late the temperature, we use the local ionized fraction as a 
proxy for the temperature (which has a few uncertainties). 
At very high redshifts things are reasonably simple since it 
is typically a simple transition from largely neutral to al- 
most completely ionized as the ionization front overruns the 
cell in question. A noted problem is what to do after this 
point. Since the ionization suppresses the low-mass sources, 
some regions start to recombine (and cool), which eventu- 
ally should allow the formation of low-mass sources again, 
but the question is when does this occur? This is a source of 
uncertainty and so a number of cell fraction thresholds need 
to be tested in future work. At the moment, the cell fraction 
(xthresh) was set to 0.1 for this work but arguments can be 
made to raise this threshold to higher values of Xthresh — 
0.5 or even 0.9 (see Appendix B of Iliev et al. 2011). Model 
1 should be retested using these thresholds to determine 
how much stochasticity one would expect from varying this 
threshold. 

• Clumping; volume averaged recombination rates in an 
inhomogeneous IGM scale with the clumping factor C = 
(p 2 ) I (p) 2 ? where the p denotes volume averaged densities. 
In our work, we set C = 1 (no clumping). The dependence 
of halo suppression for varying degrees of clumping is as yet 
unknown but will be implemented in subsequent work. 



Recently, Spitler et al. (2012) combined observations 



with novel modelling to detect evidence of inhomogeneous 
reionization by MPGCs on cluster scales. They exploited a 
fundamental characteristic of galaxy assembly (i.e. spatial 
biasing and kinematics of MPGCs) to constraint the local 
reionization epoch around individual galaxies. They found 
a joint constraint from three galaxies of z re ion = 10.5^J g 
which agrees well with the latest WMAP constraint on 
Zreion- Interestingly they also found a 1.7a indication that 
low-density environments were reionized before medium and 
high-density environments. These results are consistent with 
the theory that reionization was prolonged in duration with 
neutral-gas surviving in high-density environments for ex- 
tended periods. The redshift range of reionization of the 
local IGM in this work is consistent with their results to 
within la. 

Only recently, have computational developments made 
possible the ability to simulate vast cosmological volumes 
with exquisite resolution. Future studies such as the work in 
this paper could be applied to high resolution cluster envi- 



ronments to verify the claims made by Spitler et al. (2012) 



Preliminary work will be carried out on the Millennium-II 
simulation but more tailored simulations for this type of 
study would be desirable. 
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6 CONCLUSIONS 

Over the course of this work, we have drawn a number of im- 
portant conclusions regarding the reionization of the inter- 
galactic medium by metal-poor globular clusters. These con- 
clusions are as follows: 

(i) Treating the ionization field in an inhomogeneous 
manner results in substantially different numbers and spatial 
distributions. Our basic model adopting this type of treat- 
ment produces comparable numbers of GC candidates to 
that of the Milky Way but results in radial distributions 
which are too shallow. Whilst this model does not rule out 
possibility of GCs forming inside dark matter halos at high 
redshift, it does indicate that GCs must have formed via 
other mechanisms to account for the lack of numbers at low 
radii in our model. 

In our so-called 'delayed model' where GCs are able to 
become active after having been previously suppressed, the 
number of candidate GCs drastically increases. These addi- 
tional candidates however do not end up residing within the 
inner, R < 60 kpc, of the host halo at z = and so the re- 
sulting radial distributions are still not comparable to that 
of the Milky Way. 

In the arbitrary truncation model, the radial profiles and 
numbers are comparable to the Milky Way metal-poor GCs. 
This is because in the truncation model (a) there are extra 
sources available at high redshift which boosts the number 
of sources at small radii and (b) it by definition, totally 
removes the sources below ztrunc which reduces the number 
of sources at large radii. This ultimately shows that realistic 
populations can be produced if external sources dominate 
the photon budget at later times. 

Overall, the more sophisticated work in this paper clearly 
shows that treating the ionization of the local IGM in a 
spatially inhomogeneous manner leads to the IGM being 
reionized at different places at different times which greatly 
impacts the radial distributions of MPGCs. 

(ii) Globular clusters injected high numbers of ionizing 
photons into the IGM at high redshift. Conservative effi- 
ciency estimates of f 1 = 150 photons/baryon would have 
resulted in MPGCs ionizing more than 50% of the local IGM 
mass and 60% of the local IGM volume by z = 10 (within 
an enclosed volume of 2 3 h~ 3 Mpc 3 centred on the host) 
and upwards of 98% of the volume and 90% of the mass 
for the photon-richest model (/ 7 = 5000 photons/baryon). 
We also estimate that as many as 10 70 photons (minimum: 
10 68 - 5 photons) were injected in the IGM from MPGCs alone 
during the early build up the Milky Way. Such a quantity 
of photons would surely impact the formation of structures 
around not only the Milky Way, but other GC rich environ- 
ments as well. 

(iii) The suppression rate of MPGCs in our simplest 
reionization model (Ml_512_phl50) was 36 % of a total 
310 possible MPGC sources. The number of objects which 
form between the photon-poor (/ 7 = 150) and photon-rich 
(/ 7 = 5000) environments differs by a factor of 2.7. 

(iv) The unsuppressed MPGCs in all models have a nar- 
row age range (mean = 13.34 Gyr, a = 0.04 Gyr) consistent 
with current ages estimates of the Milky Way MPGCs. 

(v) The radial distributions for the Ml_512_phl50 and 
Ml_512_ph250 models were the only distributions to be sta- 



tistically consistent with the distribution of the Milky Way 
MPGCs at the 1% significance level. 

(vi) In an extension to Model 1, we found that dynamical 
destruction destroys nearly half (52%) of the total GCs avail- 
able resulting in a sample of GC candidates which is much 
smaller in number when compared to the Milky Way metal- 
poor GCs. A KS-test carried out on all six Model 2 distribu- 
tions rejected the null hypothesis that they are drawn from 
the same distribution as the Milky Way MPGCs at the 1% 
significance level. 

(vii) Allowing suppressed halos to turn on once they sat- 
isfy the required environmental conditions results in double 
and some times triple the number of potential GCs having 
formed. These 'delayed GCs' begin forming at z ~ 19 and 
continue to form up until the present day. The majority of 
these objects however form within neutral-gas regions of the 
simulation volume at z < 8 and at R < 300 kpc from 
the host. Due to their late formation times, these objects 
are most likely self-enriched and are therefore not suitable 
for comparison to the Milky Way MPGCs. If we impose 
truncation at z = 10, the photon-rich radial distributions 
(/ 7 = 1000 — 5000 photons/baryon) are consistent with that 
of the Milky Way. 

In summary, a number of avenues could be explored in 
future work. Measurements of the variation in our results 
across a wide variety of Milky Way-type halos and a number 
of other sources of reionization should be studied. Further- 
more, testing a range of formation efficiencies, cell fraction 
thresholds (xthresh) and an improved treatment of the dy- 
namical disruption will help comprehend the ultimate origin 
of GCs and the role they played during the epoch of reion- 
ization. 
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